function [RX, UA] = DSP_BPS_SP(RX, varargin)
% parallel BPS

    p = inputParser;
    p.addRequired('M');
    p.addParameter('NTestPhases', 128);
    p.addParameter('BlockSize', 20);%depending on the strength of phase noise
    p.addParameter('NWorkers', 1, @(x) x>0 && floor(x)==ceil(x));
    p.parse(varargin{:});

    M   = p.Results.M;
    BS  = p.Results.BlockSize;
    NTP = p.Results.NTestPhases;
    NW  = p.Results.NWorkers;

    RX = RX(:);

    % normalize signal
    RX = RX/sqrt(mean(abs(RX).^2))*sqrt(2*(M-1)/3);
    % RX = RX/sqrt(mean(abs(RX).^2));

    if NW == 1 % serial processing
        A = dispatch(RX, M, NTP, BS);
    else       % parallel processing
        fprintf('parallel process, %d workers are used\n', NW);

        L      = numel(RX);
        LParts = fix(L/NW);
        OL     = fix(BS/2);
        RX_    = {};

        parfor n = 1:NW
            % add overlaps
            if n < NW, end_ = n*LParts+OL; else, end_ = L; end
            if n > 1, start_ = (n-1)*LParts+1-OL; else, start_ = (n-1)*LParts+1; end
            RX_{n} = RX(start_:end_);

            % dispatch jobs
            A_{n} = dispatch(RX_{n}, M, NTP, BS);

            % remove the overlaps
            if n < NW, A_{n} = A_{n}(1:end-OL); end
            if n >  1, A_{n} = A_{n}(OL+1:end); end
        end

        A = reshape([A_{:}], [], 1);
    end

    % unwrap phase
    UA = unwrap(4*A)/4;
    RX = RX.*exp(1j*UA);
    RX=RX.';
end

function A = dispatch(RX, M, NTP, BS)
    L = numel(RX);
    NBlocks = L-BS+1;

    RX = reshape(RX, 1, []);

    P = zeros(1, L);
    A = zeros(1, L);

    testAngles = (0:NTP-1)*pi/(2*NTP);
    testPhases = exp(1j*testAngles).';

    % initial block evaludation
    for i = 1:1
        mid      = i+fix(BS/2);
        W_i      = RX(i:i+BS-1);
        tW       = testPhases * W_i;
        dW       = DSP_DecisionQAM_SP(tW, M);
        err      = sum(abs(dW - tW).^2, 2);
        [~, idx] = min(err);
        P(mid)   = testPhases(idx);
        A(mid)   = testAngles(idx);
    end

    % slided processing
    for i = 2:NBlocks
        mid      = i+fix(BS/2);
        C_i      = RX(i+BS-1);
        tW       = [tW(:,2:end), testPhases*C_i];
        dW       = [dW(:,2:end), DSP_DecisionQAM_SP(tW(:,end), M)];
        err      = err - (abs(dW(:,1) - tW(:,1)).^2) + (abs(dW(:,end) - tW(:,end)).^2);
        [~, idx] = min(err);
        P(mid)   = testPhases(idx);
        A(mid)   = testAngles(idx);
    end

    % populate the leading and trailing phase entries
    P(1:fix(BS/2)-1)         = P(fix(BS/2));
    P(NBlocks+fix(BS/2)+1:L) = P(NBlocks+fix(BS/2));
    A(1:fix(BS/2)-1)         = A(fix(BS/2));
    A(NBlocks+fix(BS/2)+1:L) = A(NBlocks+fix(BS/2));

    % compensate the input signal
    A  = A(:);
end
